Tension and stiffness of the hard sphere crystal-fluid interface 
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A combination of fundamental measure density functional theory and Monte Carlo computer simulation 
is used to determine the orientation-resolved interfacial tension and stiffness for the equilibrium hard-sphere 
crystal-fluid interface. Microscopic density functional theory is in quantitative agreement with simulations 
and predicts a tension of 0.66 ksT/a^ with a small anisotropy of about 0.025 fceT and stiffnesses with 
e.g. 0.53kBT/a^ for the (001) orientation and l.OSfcBT/a^ for the (111) orientation. Here fceT is denot- 
ing the thermal energy and a the hard sphere diameter We compare our results with existing experimental 
findings. 
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Solidification and melting processes involve crystal-fluid 
interfaces that separate the disordered from the ordered phase. 
Understanding the properties of such interfaces on a micro- 
scopic scale is pivotal to control and optimize crystal nucle- 
ation and the emerging microstructure of the material. Impor- 
tant applications include the fabrication of defect-free metal- 
lic alloys Uj] and of photonic |2(], phononic yy and protein 
141] crystals. In equilibrium, i.e. between a coexisting crystal 
and fluid phase, creating a crystal-fluid interface results in a 
free energy penalty per area that is called interfacial tension. 
Unlike the liquid-gas or fluid-fluid interface, the structure of 
the solid-fluid interface depends on its orientation 151]. This 
anisotropy is associated with a difference between the inter- 
facial tension and the interfacial stiffness of a crystalline sur- 
face. 

Predicting crystal-fluid interfacial tensions by a molecular 
theory is a very challenging task. Classical density functional 
theory of freezing provides a unifying framework to describe 
the solid and liquid on the same footing and is therefore in 
principle a promising tool. In this respect, the simple ather- 
mal hard sphere system which exhibits a freezing transition 
from a fluid into a face-centered-cubic (fee) crystal, is an im- 
portant reference system. The accuracy of previous density 
functional calculations of the hard sphere solid-fluid interface 
IQ^, however, was hampered by the lack of knowledge of a 
reliable functional and severe restrictions in the parametriza- 
tion of the trial density profile. 

In this letter, interfacial tensions and stiffnesses of the equi- 
librium hard sphere crystal-fluid interface are predicted us- 
ing fundamental measure density functional theory lIlOll which 
has been shown to predict accurate bulk freezing data Ollll . 
The interfacial tension and stiffness for five different orien- 
tations are obtained, namely along the (001), (Oil), (111), 
(012) and (112) orientations (see Fig. ^. A small orienta- 
tional anisotropy for the tensions is found and the average 
tension is about 0.66 k-^T /a"^ with k^T denoting the ther- 
mal energy and a the hard sphere diameter. For the stiff- 
nesses the data are spread in a much wider range between 



0.28/08^/(7^ for the (Oil) orientation with lateral direction 
[100] and I.OS/sbT/ct^ for the (111) orientation. We have 
also conducted Monte Carlo simulations to extract the stiff- 
ness from capillary wave fluctuations for the above orienta- 
tions except (012), thereby improving the accuracy of earlier 
data II12I - II8II . We find quantitative agreement between density 
functional theory and computer simulation. 

In equilibrium, the athermal hard sphere model system is 
solely characterized by the volume fraction 0; the thermal 
energy k-oT just sets the energy scale. The fluid-solid(fcc) 
freezing transition is first-order with coexisting fluid and solid 
volume fractions of 0f = 0.492 and (j)s = 0.545, respectively, 
and a coexistence pressure of Pc = 11-576 k^T /a^ i\% . For 
a given volume V containing coexisting bulk fluid and sohd, 
and a planar fluid-solid interface of area A, the excess grand 
free energy per area is the surface or interface tension, given 
by 7 = (51 + pV)/A, with fi denoting the grand-canonical 
free energy. For crystal-fluid interfaces, 7 depends on the ori- 
entation of the interface, characterized by a normal unit vector 
n relative to the crystal lattice. The latter is fixed with the fee 
cubic unit cell edges parallel to the Cartesian coordinate axes 
of our system, see Fig.[T] 

The central quantity to describe thermal fluctuations, 
i.e. capillary waves, along the anisotropic crystal-fluid inter- 
face is the interfacial stiffness defined tensorially 112011 as 



lafiin) ="f{n) 



dfiadfiis 



(1) 



for two directions fia and hfi that are orthogonal to h. 

We calculate the tension of the hard sphere crystal-fluid in- 
terface using classical density functional theory (DFT) that 
prov ides direct access to the grand-canonical free energy Q, 
1 2111 . We employ the geometric fundamental measure ap- 
proach first established by Rosenfeld 11221 12311 and most ac- 
curately refined in the so-called White Bear version mark II 
. A free minimization of this theoiy in the bulk phases 
gives accurate hard sphere bulk coexistence data which 
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are needed as a reliable input for the calculation of interfacial 




FIG. 1 : (Color online) In the left panel, the surface orientations, as 
listed in Table U are indicated on an octant of the unit sphere. The 
right panel shows a Wulff plot of the corresponding interfacial ten- 
sion 7(n) ; here, the colors display the value of the tension for a given 
orientation. 
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tensions. The crystal-fluid phase transition occurs at a coex- 
istence chemical potential jic/k-QT = 16.3787 and a coexis- 
tence pressure pccr^ /k-oT ~ 11.8676. The coexistence pack- 
ing fractions of the fluid and solid are respectively (j)f = 0.495 
and (t>s = 0.544, in close agreement with the aforementioned 
computer simulation data 11911 . 

At the prescribed coexistence chemical potential /ic, the 
grand free energy functional is numerically minimized inside 
a rectangular cuboid box of lengths L^., Ly, and L^ with pe- 
riodic boundary conditions in all three directions |J7|. The 
surface normal is pointing along the z-direction and the box 
length Lz is chosen large enough (about 50 — 60a) to ensure a 
large part of bulk crystal and fluid phase at coexistence which 
are separated by two interfaces. The lateral dimensions Lx 
and Ly of the box depend on the surface orientation relative 
to the fee crystal. They are determined by the minimal size 
of a periodic rectangular cross section which accomodates the 
prescribed relative orientation. The density field is resolved 
on a fine rectangular grid in real space with a spacing of about 
0.02(7. Starting from an initial profile which contains the 
two bulk parts of pre-minimized crystal and fluid, the density 
functional is minimized using a Picard iteration scheme com- 
bined with a direct inversion in the iterative subspace method 
|24|,|25|] and a simulated annealing technique ^M- Finite size 
effects due to the finite grid size were excluded by also us- 
ing smaller grid spacings to ensure free minimization of the 
density functional in practice. 

Results for the minimized density profiles are dis- 
played in Fig. |2] for five different orientations. Both 
the laterally integrated (z-resolved) density field p{z) = 

L £ /o "^ /o " P^^' y^ z)dydx and contour plots, p{x = 
0, y, z), are shown. 

The DFT results for the interfacial tension are summarized 
in Table U for five different orientations. With a slight orien- 
tational dependence, all the values vary around 0.66 k^T/a^. 
The errors given in Table I] are estimated from several inde- 
pendent minimization runs. Since the anisotropy is weak, the 
orientational resolved interfacial tension can be well-fitted by 
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FIG. 2: (Color online) DFT results: a) Laterally integrated density 
profiles p{z) for the five surface orientations, as indicated; b) contour 
plots at a:: = 0. The periodic length of the total profiles in z-direction 
is 50.15a (001), 53.19a (Oil), 65.15a (111), 56.07a (012), and 
61.42a (112). 
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with h = (ni, 71,2, na), Q = rif + n| + n^, S 
four fit parameters 70, ei, 62, £3. The expansion ^ can be 
used to obtain the interfacial stiffness {B from the DFT data 
of the anisotropic interfacial tension 11 19[ 12711 . The resulting 
anisotropy of the stiffness is considerably larger than the one 
of the tension (thus, the resulting errors of 7^^ (h) are also 
larger). For the five different orientations considered in this 
work, the resulting data for the interfacial stiffness and for the 
fit parameters are listed in Table H] 

In the Monte Carlo (MC) simulations, similar to the pro- 
cedure in II12I [191], inhomogeneous hard-sphere systems at 



the coexistence pressure pc are prepared followed by pro- 
duction runs in the canonical ensemble. The canonical MC 
simulation consists of particle displacement moves according 
to a standard Metropolis criterion where the trial displace- 
ments of the particles are randomly chosen from the interval 
[-0. 1 1 a,+Q. 1 Icr]. The inhomogeneous solid-fluid systems are 
placed in rectangular cuboid simulation boxes of nominal size 
i X L X 5L (L « 25 a), containing about 10'^ particles. We ap- 
ply periodic boundary conditions in all three dimensions, the 
fee crystal with z-extension of about 2L resides in the mid- 
dle of the box and is separated from the fluid by two planar 
interfaces. Since the system is in equilibrium, the amount of 
crystal and fluid phase as well as the interfaces remain stable 
during the simulation. We consider the crystal orientations 
(001), (01 1), (1 11) and (1 12), see Fig.[I] At each orientation, 

10 independent runs are performed and in each of these runs, 
10,000 independent configurations are generated that are used 
for the analysis of the interfaces. 

The stiffnesses 7 are extracted from the capillary wave 
spectrum {h'^(q)) 112711 . i.e. the correlation function of the in- 
terface height fluctuation h{q) (with q = (q^^qy) the two- 
dimensional wave-vector along the lateral extension of the in- 
terface). In order to determine h(q) a criterion has to be intro- 
duced according to which one can di stinguis h between fluid 
and crystal particles. Following work lll2[ll9ll . the rotational- 
invariant bond-order parameter qeqeii) was used 128112911 . To 
distinguish between crystalline and fluid particles, we adopt 
the same criterion as in Refs. lll2[ll9ll . where a particle i was 
identified as one with crystalline order if qeqeii) > 0.68, oth- 
erwise it was defined as a liquid-like particle. Moreover, the 
local position of the interface is defined by the set of crys- 
talline particles at the interface (particles which have less than 

1 1 crystalline neighbors). Some particles in the liquid bulk 
identified as crystalline were removed by searching the largest 
cluster among the particles identified as interface-particles. 
The fluctuation of the local interface position is defined as 
h{xi,yi) — Zi — (z), with i the index of a particle on the 
surface and (z) the instantaneous average location of the inter- 
face. The irregularly distributed particle coordinates (xi, yi) 
are then mapped onto a regular grid in the xy plane with grid 
spacing Ax = Ay = a using Shepard interpolation ll27ll . Fi- 
nally, the height fluctuation ft.(g) is obtained from a Fourier 
transformation of the interpolated heights. 

Figure [5] shows the g-dependent interfacial stiffness, 
as defined by the equation ^i(qx)ql + J2{qy)ql = 
kBT/[L^Ly{h^{q))]: for the (001) and (111) orientation 
7(9) = 71(92^) = l2{qy) holds whereas for the (Oil) and 
(112) orientation there are two different coefficients ji(qx) 
and 72 (^Zy) that can be determined from the latter equation by 
considering qy = or q^ = 0, respectively. The solid lines in 
Fig. [3] are fits of the data for q < 1.5 a^^ with the function 
7(9) = 7 + o,q^ + bq* yielding the values for the stiffness 7 
for g — > 0. 

In TableHl the values of 7, as obtained from our simulation, 
are given in comparison to previous simulation results and to 
DFT. A direct comparison is not possible since in DFT the 
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FIG. 3: (Color online) g-dependent interfacial stiffness 7(g) for the 
(001) and (ill) orientation [a)] as well as the (Oil) and (112) orien- 
tation for the indicated directions [b)]. Note that for (112) only the 
[110] direction is shown because 7(g) for the [111] direction is very 
similar to that of the latter direction. 



TABLE I: Interfacial tensions 7 and stiffnesses 7 in units of ksT/a^ 
for different surface normal vectors (round brackets) and tangential 
directions (square brackets). In DFT the tensions are measured di- 
rectly, in the simulation the stiffnesses. The other data are listed 
italicized and are calculated using the fit function (|2j. The fit pa- 
rameters are obtained from a least square fit to the measured data. 
For DFT they are 70 = 0.664(2)fcBT/cr^ ei = 0.1076(120), 
£2 = -0.01364(292), £3 = -0.0023(209) and for simulation 
70 = 0.618(ll)fcB^/c^^ £1 = 0.0905(32), £2 = -0.00547(44), 
£3 = 0.0054(25). As a reference previous simulation results for 
tensions |17ll and stiffnesses 11511 are shown in the last column. The 
numbers in parentheses indicate the uncertainty in the last digit(s). 





orientation 


theory 


simulation 


Ref. |T5,iT] 


7 
7 


(001) 
(001) 


0.687(1) 
0.53(14) 


0.639(11) 
0.419(5) 


0.5820(19) 
0.44(3) 


7 
7 
7 


(Oil) 

(011)[100] 

(011)[011] 


0.665(1) 

0.283(35) 
0.86(14) 


0.616(11) 
0.401(5) 
0.769(5) 


0.5590(20) 
0.42(3) 
0.70(3) 


7 
7 


(111) 
(111) 


0.636(1) 
1.025(79) 


0.600(11) 
0.810(5) 


0.5416(31) 
0.67(4)" 


7 
7 
7 


(012) 

(012)[100] 

(012)[021] 


0.674(5) 
0.454(57) 
0.71(12) 


0.623(11) 
0.575(5) 
0.420(5) 


0.5669(20) 
0.59(3) 
0.43(3) 


7 
7 
7 


(112) 

(112)[110] 

(112)[111] 


0.654(1) 
0.973(41) 
0.704(50) 


0.611(11) 
0.606(5) 
0.550(5) 





" This value is for the rhcp-crystal-liquid, rather than the 
fcc-crystal-liquid interface. See U^l for details. 



tensions are calculated whereas in MC the stiffnesses are mea- 
sured. A comparison is only possible using a tension-stiffness 
conversion through a least-square fit to the tension anisotropy 
expansion ^ and the coiTesponding expression for the stiff- 
nesses (through a combination of (|2]i and ([T]i), giving the av- 
erage tension 70 and the parameters ti (z=l,2,3). [Because the 
fit function ^ can not reproduce the inner anisotropy for the 
(012) orientation (shown in the theory column of Tab. H) we 
have not taken into account the simulation data for the stiff- 
nesses at the orientations (012) and (112) for the least-square 



fit. The different inner anisotropy therefore is not a shortcom- 
ing of DFT.] Here, an element of uncertainty is added by the 
truncation of the expansion since the single terms especially in 
the stiffness expansion are not small (note also the associated 
error bars in converted quantities). 

As expected for a fcc-fluid interface, DFT shows the largest 
interfacial tension for the (001) interface orientation and 
the lowest one for the (111) orientation, giving the tension 
anisotropy (7(001) - 7(lll))/2 = 0.025fcBT/ cr^. The av- 
erage tension 70 = 0.664 /cbT/ct^ is 7.4% higher than that 
from the simulation. This deviation most likely stems from 
the fact that in DFT (long-ranged) fluctuations in the interface 
are averaged out. A comparison between the stiffnesses shows 
deviations from up to 0.36 for the (112) [110] direction to less 
than one percent for the (012) [021] direction. 

Previous simulations obtained the values 
0.559(17) fceT/cr^ |[ll and 0.5610(12) fceT/cr^ [[13] for 
7o. These values are 10% smaller than our simulation results. 
An obvious discrepancy appears for the (111) interface 
orientation where we have not observed deviations from a fee 
packing in contrast to fll5ll . Further differences to previous 



simulations are the use of a different geometry and of a 
rotational-invariant order parameter for the identification of 
crystalline particles. 

We now compare our data to real-space experiments on dis- 
persions of hard-sphere-like colloids. They often carry resid- 
ual charges and are polydisperse. This renders a comparison 
with theoretical results on hard spheres difficult. Hitherto the 
interfacial tension was indirectly measured by interpreting the 
probability to find small (non-spherical) clusters in terms of 
classical nucleation theory OOll yielding data for a mean ten- 
sion of about 7o =0.1 k-oT/a^. Clearly, given the limitations 
of the hard sphere model due to particle charging and the in- 
herent assumption of small spherical crystalline nuclei, this is 
just a rough estimate of 79. An alternative experimental route 
is via the analysis of the capillary-wave spectrum similar to 
what we do in our MC simulations 13 11 - 13311 pro viding direct 
access to the interfacial stiffnesses. In Ref. 13211 . the reported 
stiffness of 1.2k'QT/a^ for an interface between a randomly 
stacked hexagonal close packed (hep) crystal and its melt is 
significantly higher compared to our results which might re- 
flect the slight charge, the limited ensemble averaging, and an 
ad-hoc value for the viscosity required for the analysis in this 
experiment. In Ref. 113 IL the reported stiffnesses were in the 
range of about (0.7 - l.'i)k^T /a'^. Interestingly, the stiff- 
ness for the (Oil) interface was found to be isotropic and the 
highest value for the stiffness was found for the (001) orien- 
tation, oppositely to what the authors expected and what we 
found for hard spheres. This might be due to a limited num- 
ber of crystalline layers and the small gravitational length of 
ct/7 which limits the thickness of the liquid. Finally, Nguyen 
et al. i33ll grew crystals of PNIPAM particles in a tempera- 
ture gradient and analyzed the capillary waves along crystal- 
fluid interfaces after the removal of the temperature gradient. 
They measured averaged stiffnesses for several interface ori- 
entations in the range of (0.19 — 1.13)fcBT/(T^ that show the 



best agreement with our results. Nevertheless, the latter ex- 
periment is also not accurate enough to validate theory and 
simulation on a quantitative level and thus more experimental 
explorations are required. 

In conclusion, we have predicted accurate values for the 
anisotropic crystal-fluid surface tensions and stiffnesses of a 
hard sphere system by using both fundamental measure den- 
sity functional theory and Monte Carlo simulations. A small 
anisotropy in the tensions of about 10% was found which is, 
however, crucial for the transformation to stiffnesses which 
differ up to a factor of 4. These predictions can help to clar- 
ify apparent discrepancies found in real-space experiments of 
sterically-stabilized colloidal suspensions 130143 3ll . Since the 
anisotropic tensions control changes of the interfacial shape, 
their precise quantitative determination help to understand 
crystal nucleation 11341 13511 and the transport of larger carri- 
ers through the interface. They may also serve as further input 
to phase-field-crystal calculations which explore solidification 
processes on larger length and time scales 113 61 - 13811 . 

Future work should address soft interactions and attractions 
(as relevant, e.g. for colloid-polymer mixtures), in order to 
scan the full range from a fluid-crystal to a vapor-crystal in- 
terface. Further extensions can be done along similar ideas 
as used and proposed here for binary mixtures. Finally the 
recent extension of DFT towards dynamics for Brownian sys- 
tems can be used to explore the time-dependent growth ki- 
netics and relaxation towards equilibrium for the hard sphere 
interface ll39ll . 
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